/*
 * main.cpp
 *
 *  Created on: 26.02.2013
 *      Author: Roma Jam
 */
#include "time.h"
#include "fft_rj.h"
#include "main.h"
#include "complex"

using namespace std;

#define N0          9           // dlina registra
#define K           (1 << N0)   // dlina PSP
#define SIZE        K         // Period of search

#define DopCount    5000

int aPSP[SIZE]; // PSP
unsigned int RG = 0xffffffff;
unsigned int O5 = (1 << 5);
unsigned int O9 = (1 << 9);
unsigned int O10 = (1 << 10);
complex<double> cPSP1[SIZE],
                cPSP2[SIZE],
                cPSP3[SIZE],
                cPSP4[SIZE];
fftw_plan Dft_of_PSP;
double Fd_min = 0;
double Fd_max = 1000;

double Fdop[DopCount];
double Fd;
double Fref = 0;
double dT = 2e-6;
double Ttek = 0;
double r1 = 0;
double r2 = 0;
double Amp = 0.5;

int FreqNow;
int MaxIndex;
FILE *Txt, *DopArr;

int main(void){
    Txt = fopen("PSP.txt", "wt");
    srand((int)time(0));
    printf("Start\n");


    for (FreqNow = 0; FreqNow < DopCount; ++FreqNow) {
        Fd = FreqNow;
        GeneratePSP();

        Transform.FFT(cPSP1, cPSP1, SIZE, tf_forward);
        Transform.FFT(cPSP2, cPSP2, SIZE, tf_forward);

        for (int i = 0; i < SIZE - 1; i++) {
            cPSP1[i] = conj(cPSP1[i]);
            cPSP2[i] *= cPSP1[i];
        }

        Transform.FFT(cPSP2, cPSP2, SIZE, tf_backward);
        Transform.Norm(cPSP2, SIZE);

        MaxIndex = Transform.MaxSearch(cPSP2, SIZE);
        Fdop[FreqNow] = real(cPSP2[MaxIndex]);

        if(FreqNow == 1){
            printf("Saving PSP\n");
            for (int i = 0; i < SIZE; ++i) {
                fprintf(Txt, "%i %lf %lf\n",
                        aPSP[i],
                        real(cPSP2[i]), imag(cPSP2[i])
                );
//                printf("%lf + j%lf\n", real(cPSP2[i]), imag(cPSP2[i]));
            }
        }
        Transform.Reset(cPSP1, SIZE);
        Transform.Reset(cPSP2, SIZE);
    }

    DopArr = fopen("Dop.txt", "wt");
    for (int f = 0; f < DopCount; ++f) {
        fprintf(DopArr, "%lf\n", Fdop[f]);
    }
    fclose(DopArr);

    printf("Ok!\nFile had been saved\n");
    fclose(Txt);
}

void GeneratePSP(){
    int k;

    Ttek = (double)rand() / RAND_MAX * 12357913.0;
    ModPSP();

    for(k = 0; k<SIZE-1; k++){
        r1 = Random();
        r2 = Random();

        cPSP1[k] = (double)aPSP[(k+100)%K] * Amp * exp(complex<double>(0.0,-2.0*M_PI*Fd*Ttek));
//                 + (double)aPSP[(k+400)%K] * exp(complex<double>(0.0,-2.0*M_PI*Fd*Ttek));
//                 + (complex<double>(1*sqrt(-2.0*log(r1))*cos(2.0*M_PI*r2),1*sqrt(-2.0*log(r1))*sin(2.0*M_PI*r2)));

        cPSP2[k] = (double)aPSP[k] * exp(complex<double>(0.0,-2.0*M_PI*Fref*Ttek));   // form PSP2

        cPSP3[k] = (double)aPSP[k] * exp(complex<double>(0.0,-2.0*M_PI*Fref*Ttek));         // real PSP3


        Ttek += dT;
    }
}

void ModPSP(){
    for(int k = 0; k<SIZE-1; k++){
        RG = (RG & 0xfffffffe) | (((RG & O5)>>5)^((RG & O9)>>9));
        RG <<= 1;
        aPSP[k] = (RG & O10)>>10;
        if(aPSP[k] == 0) aPSP[k] = -1;
    }
}

double Random(){
    double r;
    do{
        r = (double)rand() / RAND_MAX ;
    } while(r == 0.0);
    return r;
}

// Nikolay Klimentkin
